#012 - Portfolio Optimization II - Log Returns¶

In this code we'll optimize code used in #011, try to peform a greater number of simulations using log returns in a continuous function, and then convert the results only to arithimetical returns. Expected to reduce considerably time spending in monte carlo runs.

simplified calcs using risk free = 0, makes sharpe = rm/std

libs¶

In [1]:
#import Libraries
import pandas as pd
from pandas_datareader import data as pdr
import numpy as np
from scipy.optimize import minimize
import random

import plotly.express as px
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick

Functions¶

In [2]:
def rand_weights(n):
    ''' Produces n random weights that sum to 1 '''
    k = np.random.rand(n)
    return k / sum(k)
In [3]:
def run_portfolio(df, weights, initial_investment, risk_free):
    ''' Performs an asset allocation and calc Portfolio Metrics'''
    # The function returns: 
    # (0) Portfolio Data Frame 
    # (1) Expected portfolio return 
    # (2) Expected volatility 
    # (3) Sharpe ratio 
    # (4) Return on investment 
    # (5) Final portfolio value in dollars
    
    # Perform asset allocation using the random weights (sent as arguments to the function)
    portfolio_df =  df.copy()
    # Scale stock prices using the "price_scaling"
    scaled_df = df.copy()
    for i in df.columns[0:]:
        scaled_df[i] = df[i]/df[i][0]
    #enumerate method links Stocks tickers in columns along with a counter position weight (i), like an index
    for i, stock in enumerate(scaled_df):
        portfolio_df[stock] = weights[i] * scaled_df[stock]  * initial_investment
    
    # Sum up all values and place the result in a new column titled "portfolio value [$]" 
    portfolio_df['Total Value [$]'] = portfolio_df.sum(axis = 1, numeric_only = True)
    # Calculate the portfolio percentage daily return and replace NaNs with zeros
    portfolio_df['Daily Return [%]'] = portfolio_df['Total Value [$]'].pct_change(1) * 100 
    portfolio_df.replace(np.nan, 0, inplace = True)
    # Calculate the return on the investment = last final value of the portfolio compared to its initial value
    return_on_investment = ((portfolio_df['Total Value [$]'][-1:] - portfolio_df['Total Value [$]'][0]) 
                            / portfolio_df['Total Value [$]'][0]) * 100
    
    # Daily change of every stock in the portfolio 
    portfolio_daily_return_df = portfolio_df.drop(columns = ['Total Value [$]', 'Daily Return [%]'])
    portfolio_daily_return_df = portfolio_daily_return_df.pct_change(1)
    portfolio_daily_return_df.replace(np.nan, 0, inplace = True)  
  
    # Portfolio Expected Return formula
    expected_portfolio_return = np.sum(weights * portfolio_daily_return_df.mean() ) * 252
  
    # Portfolio standard deviation, volatility (risk)
    covariance = portfolio_daily_return_df.cov() * 252 
    expected_volatility = np.sqrt(np.dot(weights.T, np.dot(covariance, weights)))
    
    # Calculate Sharpe ratio
    sharpe_ratio = (expected_portfolio_return - risk_free)/expected_volatility 
    
    return expected_portfolio_return, expected_volatility, sharpe_ratio, portfolio_df['Total Value [$]'][-1:].values[0], return_on_investment.values[0]

11.1 Read and Organize Data and Run Simulations¶

In [4]:
# file_name =          input('Input the CSV file name: ')
# initial_investment = int(input('Input the initial investment: '))
# risk_free =          float(input('Input the risk free annual rate: ')) # https://ycharts.com/ind icators/10_year_treasury_rate 
# sim_runs =           int(input('Input how many runs to perform on Monte Carlo Simulation: '))

file_name =          'BRL5'   # CSV file with 5 Brazilian stocks and 10 year Ajusted Close Prices
initial_investment = 1000000    # 1mm
risk_free =          0.039      #3.9%
runs =               100000    # 100k
In [5]:
# read CSV file
df = pd.read_csv(file_name)
# replace the colum Date to Index
df.set_index(['Date'], inplace = True)

tickers = list(df.columns)

11.2 Portfolio Metrics and Run Monte Carlo Simulations¶

In [6]:
#Log is continuous and approx arithmetic daily return (252. over basis)
log_returns = df.pct_change().apply(lambda x: np.log(1+x)).dropna() #return for each stocks position column 
mean_log_returns = log_returns.mean()          # mean of Log returns | Mu
covariance_matrix = log_returns.cov()          # covariance matrix   | Sigma
In [7]:
#create arrays to store results
log_expected_return_array = np.zeros(runs)        # expected return  | E(R)
volatility_array = np.zeros(runs)                 # variances matrix | Sigma
sharpe_array = np.zeros(runs)                     #    
weights_array = np.zeros((runs, len(df.columns))) #                  | W
In [8]:
# loop to perform monte carlo with random weights
for i in range(runs):
    w = np.random.rand(len(df.columns))
    w = w/np.sum(w)
    weights_array[i, :] = w
    
    log_expected_return_array[i]     = np.sum(mean_log_returns * w * 252)
    volatility_array[i] = np.sqrt(np.dot(w.T, np.dot(covariance_matrix*252, w)))
    sharpe_array[i]              = log_expected_return_array[i]/volatility_array[i]  
In [9]:
# maximum sharpe ratio over simulation runs
opt_sharpe_index = sharpe_array.argmax()
opt_sharpe = sharpe_array.max()
opt_weights = weights_array[opt_sharpe_index]

# convert log return to arithmetic return
expected_return_array = np.exp(log_expected_return_array) - 1 

#run for optimal portfolio 
opt_return = np.exp(np.sum(mean_log_returns * opt_weights * 252)) - 1
opt_volatility = np.sqrt(np.dot(opt_weights.T, np.dot(covariance_matrix*252, opt_weights)))
In [10]:
# for y axis, select the maximum and minumum return points to plot efficient frontier
frontier_Y = np.linspace(expected_return_array.min(), expected_return_array.max(), 50)

# select each allocation, calculate and covert log to arithmetic return   
def select_returns(w_1): 
    return_w1 = np.exp(np.sum(mean_log_returns * np.array(w_1)) * 252) - 1
    return return_w1
# check if all capital is allocated into the assets, sum shoud be equal to 1
def check_allocation(w_1): 
    return np.sum(w_1)-1
# select each allocation, calculate volatility with w_1 and covariance matrix 
def select_volatility(w_1):
    w_1 = np.array(w_1)
    w_1_volatility = np.sqrt(np.dot(w_1.T, np.dot(covariance_matrix*252, w_1)))
    return w_1_volatility

# supose a equal weighted vector to start plot
w_0 = [1/len(tickers)] * len(tickers) 
# impose bounds to plot function
bounds = tuple([(0, 1) for stock in tickers])
# set a vector to plot X axis
frontier_X = []

#run for max e min returns to plot X axis. verify if w sums to 1, if return selected appox to monte carlo return R
for R in frontier_Y: 
    # impose constraints to run least squares
    constraints = ({'type':'eq', 'fun':check_allocation}, {'type':'eq', 'fun': lambda w: select_returns(w) - R})
    # Minimize a scalar function of one or more variables using Sequential Least Squares Programming (SLSQP).
    result = minimize(select_volatility,w_0,method='SLSQP', bounds = bounds, constraints = constraints)
    frontier_X.append(result['fun'])
    
In [11]:
# Create a DataFrame that contains volatility, return, and Sharpe ratio for all simualation runs
sim_out_df = pd.DataFrame({'Volatility': volatility_array.tolist(), 'Portfolio_Return': expected_return_array.tolist(), 'Sharpe_Ratio': sharpe_array.tolist() })
In [12]:
# Plotting data
fig1 = px.line(df, x= frontier_X, y= frontier_Y)
fig = px.scatter(sim_out_df, x = 'Volatility', y = 'Portfolio_Return', color = 'Sharpe_Ratio', size = 'Sharpe_Ratio', hover_data = ['Sharpe_Ratio'], size_max=10 )
fig.add_trace(go.Scatter(x = [opt_volatility], y = [opt_return], mode = 'markers', name = 'Optimal Point', marker = dict(size=[20], color = 'green')))
fig2 = go.Figure(data=fig.data + fig1.data)
fig2.update_layout(coloraxis_colorbar = dict(y = 0.7, dtick = 3), 
                #    title = 'porfolio results ', 
                   xaxis_title = 'Volatility(STD)', 
                   yaxis_title = 'Return',
                   xaxis=dict(tickformat=".1%"),
                   yaxis=dict(tickformat=".1%"))
fig2.update_layout({'plot_bgcolor': "white"})
fig2.show()

print('Best Portfolio Metrics Based on {} Monte Carlo Simulation Runs:'.format(runs))
print("  - Portfolio Stocks                          = {}".format(tickers))
print("  - Portfolio Weights                         = {} %".format(opt_weights.round(4)*100))
print('  - Portfolio Expected Annual Return          = {:.02f} %'.format(opt_return * 100))
print('  - Portfolio Standard Deviation (Volatility) = {:.02f} %'.format(opt_volatility * 100))
print('  - Sharpe Ratio                              = {:.02f}'.format(opt_sharpe))
Best Portfolio Metrics Based on 100000 Monte Carlo Simulation Runs:
  - Portfolio Stocks                          = ['ITSA4.SA', 'PETR4.SA', 'TAEE11.SA', 'VALE3.SA', 'WEGE3.SA']
  - Portfolio Weights                         = [ 0.17  1.34 56.75  6.08 35.65] %
  - Portfolio Expected Annual Return          = 21.46 %
  - Portfolio Standard Deviation (Volatility) = 21.14 %
  - Sharpe Ratio                              = 0.92